home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_11 / jacobi1.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  1.2 KB  |  41 lines  |  [MATF/MATL]
  1. function [V,D] = jacobi1(A,epsilon,show)
  2. % [V,D] = jacobi1(A,epsilon,show)
  3. % To compute the eigenpairs of a symmetric matrix.
  4. % The original Jacobi`s method of iteration is employed.
  5. % A is an n x n matrix, input.
  6. % epsilon is the tolerance, input.
  7. % V is the diagonal matrix of eigenvectors, output.
  8. % D is the diagonal matrix of eigenvalues,  output.
  9. if nargin==2, show = 0; end
  10. D = A;
  11. [n,n] = size(A);
  12. V = eye(n);
  13. done = 0;
  14. working = 1;
  15. stat = working;
  16. cntr = 0;
  17. [m1 p] = max(abs(D-diag(diag(D)))); % Select element
  18. [m2 q] = max(m1);                   % element of largest
  19. p = p(q);                           % magnitude.
  20. while (stat==working),
  21.   t = D(p,q)/(D(q,q) - D(p,p));
  22.   c = 1/sqrt(t*t+1);
  23.   s = c*t;
  24.   R = [c s; -s c];
  25.   D([p q],:) = R'*D([p q],:);
  26.   D(:,[p q]) = D(:,[p q])*R;
  27.   V(:,[p q]) = V(:,[p q])*R;
  28.   cntr = cntr+1;
  29.   if show==1,
  30.     home; if cntr==1, clc; end; 
  31.     disp(['Jacobi iteration No. ',int2str(cntr)]),disp(''),...
  32.     disp(['Zeroed out the element  D(',num2str(p),',',num2str(q),') = ']),...
  33.        disp(D(p,q)),disp('New transformed matrix  D = '),disp(D)
  34.   end
  35.   [m1 p] = max(abs(D-diag(diag(D))));
  36.   [m2 q] = max(m1);
  37.   p = p(q);
  38.   if (abs(D(p,q))<epsilon*sqrt(sum(diag(D).^2)/n)), stat = done; end
  39. end
  40. D = diag(diag(D));
  41.